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Abstract 

X-ray  computed  microtomography  (XMT)  is  used  for  high- 
resolution,  non-destructive  imaging  and  has  been  applied 
successfully  to  geologic  media.  Despite  the  potential  of  XMT 
to  aid  in  formation  evaluation,  currently  it  is  used  mostly  as  a 
research  tool.  One  factor  preventing  more  widespread 
application  of  XMT  technology  is  limited  accessibility  to 
microtomography  beamlines.  Another  factor  is  that 
computational  tools  for  quantitative  image  analysis  have  not 
kept  pace  with  the  imaging  technology  itself. 

In  this  paper,  we  present  a  new  grain-based  algorithm  used 
for  computer  reconstruction  and  analysis  of  granular  materials 
(e.g.,  consolidated  or  unconsolidated  sands),  and  subsequent 
network  generation.  The  algorithm  differs  significantly  from 
other  methods  because  the  first  step  is  to  extract  the 
fundamental  granular  structure  from  a  3D  data  set,  which 
provides  a  wealth  of  information  such  as  grain  sizes,  aspect 
ratios,  orientations,  surface  areas,  etc.  Knowledge  of  the  basic 
granular  structure  serves  as  a  foundation  for  characterizing  the 
void  morphology  and  creating  physically  representative  pore 
networks.  The  algorithm  is  applied  to  a  sample  of  sandstone 
from  the  Frontier  Formation  in  Wyoming,  USA,  which  was 
imaged  using  Synchrotron  microtomography.  Morphologic 
and  flow-modeling  results  are  presented. 

Introduction 

Subsurface  transport  processes  such  as  oil  and  gas  production 
are  multiscale  processes.  The  pore  scale  governs  many 
physical  and  chemical  interactions  and  is  the  appropriate 
characteristic  scale  for  the  fundamental  governing  equations. 
The  continuum  scale  is  used  for  most  core  or  laboratory  scale 
measurements  (e.g.,  Darcy  velocity,  phase  saturation,  bulk 


capillary  pressure,  etc.).  The  field  scale  is  the  relevant  scale 
for  production  and  reservoir  simulation. 

Multiscale  modeling  strategies  aim  to  address  these 
complexities  by  incorporating  the  various  relevant  length 
scales.  While  pore-scale  modeling  is  an  essential  component 
of  multiscale  modeling,  quantitative  methods  are  not  as  well 
developed  as  their  continuum-scale  counterparts.  Hence,  pore- 
scale  modeling  represents  a  weak  link  in  current  multiscale 
techniques. 

The  most  fundamental  approach  for  pore-scale  modeling  is 
direct  solution  of  the  equations  of  motion  (along  with  other 
relevant  conservation  equations),  which  can  be  performed 
using  a  number  of  numerical  techniques.  The  finite  element 
method  is  the  most  general  approach  in  terms  of  the  range  of 
fluid  and  solid  mechanics  problems  that  can  be  addressed. 
Finite  difference  and  finite  volume  methods  are  more  widely 
used  in  the  computational  fluid  dynamics  community.  The 
boundary  element  method  is  very  well  suited  for  low- 
Reynolds  number  flow  of  Newtonian  fluids  (including 
multiphase  flows).  Finally,  the  lattice-Boltzmann  method  has 
been  favored  in  the  porous-media  community  because  it  easily 
adapts  to  the  complex  geometries  found  in  natural  materials. 

A  less  rigorous  approach  is  network  modeling,  which  gives 
an  approximate  solution  to  the  governing  equations.  Network 
modeling  of  a  porous  material  involves  discretization  of  the 
pore  space  into  pores  and  pore  throats.  Transport  is  modeled 
by  imposing  conservation  equations  at  the  pore  scale.  Network 
modeling  involves  two  levels  of  approximation.  The  first  is 
the  representation  of  the  complex,  continuous  void  space  as 
discrete  pores  and  throats.  The  second  is  the  approximation  to 
the  fluid  mechanics  when  solving  the  governing  equations 
within  the  networks.  The  positive  trade  off  for  these 
significant  simplifications  is  the  ability  to  model  transport 
over  orders-of-magnitude  larger  characteristic  scales. 
Consequently,  the  two  approaches  (rigorous  modeling  of  the 
conservation  equations  versus  network  modeling)  have 
complimentary  roles  in  the  overall  context  of  multiscale 
modeling.  Direct  methods  will  remain  essential  for  studying 
first-principles  behavior  and  sub-pore-scale  processes  such  as 
diffusion  boundary  layers  during  surface  reactions.  Network 
modeling  will  provide  the  best  avenue  for  capturing  larger 
characteristic  scales  (which  is  necessary  for  modeling  the  pore 
to  continuum  scale  transition). 
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This  research  addresses  one  of  the  significant  hurdles  for 
quantitative  network  modeling:  the  use  of  high-resolution 
imaging  of  real  materials  for  quantitative  flow  modeling.  We 
focus  in  particular  on  X-ray  computed  microtomography 
(XMT)  to  obtain  three-dimensional  pore-scale  images,  and 
present  a  new  technique  for  direct  mapping  of  the  XMT  data 
onto  networks  for  quantitative  modeling.  This  direct  mapping 
(in  contrast  to  the  generation  of  statistically  equivalent 
networks)  ensures  that  subtle  spatial  correlations  present  in  the 
original  material  are  retained  in  the  network  structure. 

The  specific  algorithm  described  here  was  designed  for  the 
analysis  of  unconsolidated  marine  sands,1  and  is  being  applied 
for  the  first  time  to  a  consolidated  sandstone.  The  terminology 
grain-based  algorithm  refers  to  the  fact  that  the  first  step  in 
the  algorithm  is  the  characterization  of  the  underlying  granular 
structure,  which  in  turn  is  used  to  aid  in  mapping  out  the  pore 
network  structure.  The  advantage  of  the  grain-based  approach 
is  twofold.  First,  the  identification  of  the  granular  structure 
from  a  voxel-based  image  is  less  sensitive  to  issues  such  as 
image  resolution  as  compared  to,  for  instance,  the 
identification  of  a  pore  skeleton.  Second,  knowledge  of  the 
granular  structure  provides  important  insight  into  the  locations 
and  connections  in  the  pore  space,  thus  providing  increased 
robustness  in  the  network  construction  algorithm.  A  side 
benefit  is  the  speed  of  the  algorithm. 

Background 

X-Ray  Microtomography 

X-ray  computed  microtomography  provides  nondestructive 
and  noninvasive  three-dimensional  images  of  the  interior  of 
objects  by  mapping  the  X-ray  absorption  through  the  sample. 
The  amount  of  absorption  depends  on  the  chemical 
composition  of  the  material  and  the  energy  of  X-ray.  XMT  is 
based  on  the  reconstruction  of  the  cross-section  of  an  object 
from  its  projection  data  by  passing  a  series  of  rays  through  an 
object,  and  measuring  the  attenuation  of  these  rays  using 
detectors  placed  on  the  downstream  side  of  the  object. 
Projections  are  obtained  by  measuring  the  X-ray  attenuation 
coefficients  of  the  sample  at  different  angles.  These 
attenuation  values  are  represented  in  images  as  discrete 
elements  (pixels  in  two-dimensional  images  and  voxels  in 
three-dimensional  images).  Synchrotron  radiation  has  several 
advantages  over  traditional  X-ray  sources  including  high  flux 
intensity  (number  of  photons  per  second),  high  degree  of 
collimation  (source  divergence  leads  to  image  blur),  and  the 
ability  to  tune  the  photon  energy  to  a  single  energy  or 
frequency  over  a  wide  range  using  an  appropriate 
monochromator  for  obtaining  specific-element  measurements. 

Over  the  last  two  decades,  XMT  has  played  an 
increasingly  important  role  in  the  characterization  of  porous 
media  flow  and  transport.  Due  to  its  non-destructive  nature 
and  increasingly  high-resolution,  synchrotron-based  XMT 
provides  the  high-quality  datasets  necessary  to  capture  the  3D 
microstructure  of  the  media2'4  and  the  distribution  of  fluids 
within  the  pore  space.5'7  The  ability  to  characterize  and 
correlate  the  void  space  microstructure  and  fluid  distributions 
provides  data  to  improve  and  validate  pore-scale  models.. 

Network  Modeling 


Network  modeling  has  a  long  history  in  the  oil  and  gas 
industry  beginning  with  the  landmark  paper  of  Fatt.x  For  a 
long  subsequent  period,  network  modeling  techniques 
employed  lattice-based  networks,  usually  decorated  with  a 
distribution  of  pore  (or  pore  body  in  some  literature)  and  pore 
throat  sizes.  These  lattice  based  models  are  valuable  for 
qualitative  studies  of  transport  in  interconnected, 
heterogeneous  structures.  However,  they  have  not  proved  to 
be  effective  for  quantitative  modeling  of  real  materials. 

Beginning  in  the  early  1990s,  new  techniques  were 
developed  for  quantitative  network  modeling.  Bryant  et  al.9 
created  physically  representative  network  models  from  the 
highly  characterized  Finney  packing.10  Oren  and  coworkers 
created  synthetic  (computer-generated)  sandtones  and 
developed  a  technique  for  extracting  networks  from  these 
structures.11  12  Their  group  has  continued  to  develop  this 
approach  and  the  resulting  networks  have  been  used  by  a 
number  of  other  investigators.13'15  Lindquist  et  al.  used  the 
medial  axes  of  microtomography  images  to  characterize  the 
pore  structure,16  and  this  approach  was  then  extended  to  allow 
for  direct  generation  of  network  structures.17  Thompson  and 
Fogler1*  applied  the  techniques  of  Bryant  et  al.9  to  computer¬ 
generated  packed  beds,  and  Al-Raoush  et  al.  extended  this 
work  to  allow  for  network  structures  with  arbitrary 
connectivity.19  Ioannidis  and  coworkers  have  used  simulated 
annealing  and  other  algorithms  to  produce  networks  that 
conform  to  key  statistics  in  real  materials.20'22 

For  the  network  modeling  used  in  this  work,  we  borrow 
the  terminology  physically  representative  network  models 9  to 
describe  the  general  class  of  models,  and  we  note  two 
important  characteristics  of  these  structures.  First,  these 
networks  are  a  one-to-one  mapping  of  a  specific  porous 
material  of  interest,  which  ensures  that  all  subtle  spatial 
correlations  in  the  pore  structure  are  retained.  Second,  the 
networks  are  described  using  rigorous  geometric  parameters, 
which  ensures  that  the  pore  morphology  is  not  compromised 
despite  the  need  to  discretize  the  pore  space.  This  latter  point 
contrasts  with  techniques  in  which  the  pore  structure  is 
transformed  into  a  network  of  interconnected  capillaries  from 
the  outset,  an  approach  which  has  been  shown  to  cause 
ambiguity  in  the  subsequent  modeling  of  flow.23 

It  should  be  noted  that  there  is  no  unique  or  correct 
discretization  of  most  real  pore  structures  (exceptions  being 
simples  structures  such  as  cubic  packings  of  spheres). 
Likewise,  there  is  more  than  one  approach  that  can  be  taken  to 
describe  a  network  model.  In  this  work,  we  use  the  set  of 
parameters  shown  in  Table  1  to  describe  the  network  structure. 
In  addition,  the  network  remains  linked  to  the  original 
structure  (whether  it  is  a  voxel  image  or  a  computer-generated 
material).  This  ensures  that  additional  characterization  could 
be  performed  if  warranted  by  a  particular  modeling  algorithm; 
it  would  gurantee,  in  essence,  that  none  of  the  original 
morphologic  data  are  lost. 

The  flow  modeling  itself  is  performed  by  imposing 
conservation  equations  at  each  pore  in  the  network,  which 
results  in  a  set  of  linear  or  nonlinear  algebraic  equations 
depending  on  the  physics  of  the  process  being  modeled.  A 
description  of  the  flow  modeling  algorithms  is  beyond  the 
scope  of  the  current  paper,  but  this  information  can  be  found 
in  many  other  papers  on  network  modeling. 
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■  Table  I.  Parameters  used  to  describe  the  network  structure. 


Variable 

Association 

Variable  Name 

Variable 

Type 

Dimension 

Network 

Domain  dimensions 

vector 

length 

Pore 

Location 

vector 

length 

Void  volume 

scalar 

length' 

Maximum  inscribed  radius 

scalar 

length 

Throat 

Interconnectivity  rperiodicity 

scalanvector 

Cross-sectional  area 

scalar 

length* 

Maximum  inscribed  radius 

scalar 

length 

Surface  area 

scalar 

length' 

Hydraulic  conductivity 

scalar 

length’ 

Dimension  on  hydraulic  conductivity  is  A3? 

Network  Generation  from  Voxel  Data 
Despite  the  advances  described  in  the  previous  section,  what 
remains  surprisingly  difficult  is  the  generation  of  physically 
representative  network  models  directly  from  microtomography 
images  of  real  materials.  The  difficulty  stems  from  the  distinct 
differences  in  the  form  and  scale  of  the  data  structures.  The 
XMT  data  consist  of  hundreds  of  millions  (sometimes 
billions)  of  voxels  on  a  Cartesian  grid.  In  contrast,  the  network 
representation  of  the  same  porous  medium  is  likely  to  be 
described  using  a  much  smaller  data  set  (order  1 03  —  1 05 
pores),  which  are  interconnected  via  a  complex  structure  of 
flowpaths  rather  than  aligned  with  any  specific  coordinate 
system. 

Most  previous  work  in  this  area  employs  the  medial  axis  as 
a  basis  for  characterizing  the  pore  structure.  For  discretized 
images  such  as  those  from  XMT,  the  medial  axis  is  obtained 
by  thinning  the  original  structure  until  a  one-voxel-thick 
skeleton  remains,  or  by  computing  distance  transforms  in  the 
void  space  (to  the  void/solid  interface)  and  retaining  the 
skeleton  of  local  maxima.  Different  medial-axis  structures  are 
obtained  depending  on  the  method  used,  the  order  for 
thinning,  the  type  of  distance  transform  used,  and/or  the  rules 
invoked  to  ensure  that  topology  is  retained.24  Lindquist  et  al.25 
used  the  medial  axes  of  3D  images  to  compute  statistical 
parameters  associated  with  the  void  space.  They  later  extended 
their  algorithm,  using  the  medial  axis  as  a  basis  for  obtaining 
pore  and  throat  size  distributions,  which  in  turn  were  mapped 
onto  a  network  structure.3  Specifically,  the  medial  axis  was 
trimmed  down  to  the  percolating  fraction,  and  nodes  were  then 
merged  (by  comparing  the  distance  separating  neighboring 
nodes  to  the  distance  to  local  surfaces)  to  define  pore 
locations.  Pore  throats  were  found  by  dilating  the  medial  axis 
until  the  dilated  cylinder  contacted  the  bounding  surfaces; 
pore-throat  geometries  were  then  obtained  by  triangulating 
between  the  medial  axis  and  voxels  along  the  perimeter  of  the 
constriction.  Sok  et  el.  simulated  immiscible  displacement  on 
networks  generated  by  this  algorithm  and  compared  them  to 
immiscible  displacements  on  regular  lattices  having  identical 
values  for  key  statistical  metrics.17  Mirroring  what  Bryant  et 
al.25  observed  for  single  phase  flow,  they  note  that  multiphase 
flow  behavior  is  not  reproduced  correctly  on  statistically 
equivalent  networks,  emphasizing  the  need  for  algorithms  that 
capture  true  pore  morphology.  Delerue  et  al.-0  applied  a 
medial-axis  technique  to  a  3D  image  of  a  resin-impregnated 
soil,  and  then  defined  the  pore-size  distribution  in  the  soil  by 
measuring  the  maximum  inscribed  balls  for  all  voxels 
contained  on  the  skeleton.  Mercury  intrusion  was  simulated 
directly  on  the  voxel  map  (rather  than  a  network 
representation).  Delerue  and  Perrier27  describe  in  detail  the 


various  computational  elements  used  in  the  algorithm.  Silin  et 
al.28  employed  a  similar  approach,  except  that  the  maximal 
inscribed  ball  was  found  for  each  void-phase  voxel  in  the 
packing.  Though  computationally  more  intensive,  this 
approach  allows  the  pores  to  be  found  independently  of  the 
skeleton.  They  tested  their  algorithm  using  computer¬ 
generated  sphere  packings  and  a  CT  image  of  Fontainbleau 
sandstone.  However,  pore  network  models  were  not  created, 
and  a  quantitative  assessment  of  the  pore  locations,  sizes,  and 
connectivity  was  not  made. 

Materials  and  Methods 

Sandstone  Sample 

The  sample  used  for  this  application  was  taken  from  the  Wall 
Creek  Member  of  the  Cretaceous  Frontier  Formation, 
Wyoming,  USA  as  part  of  an  integrated  geologic,  geophysical, 
and  engineering  study.29  Ten  wellbores  were  drilled  through 
the  Wall  Creek  Member  and  sample  plugs  were  selected.  The 
subject  sample  is  from  a  tidal  facies.  The  permeability  of  this 
sample  was  severely  reduced  by  calcite  cementaion.  This 
reduced  the  permeability  by  several  orders  of  magnitude  in 
approximately  15  volume  percent  of  the  tidal  facies.30  The 
sediments  preserved  in  these  rocks  were  derived  from  uplifts 
to  the  north  and  west.  Nyman  et  al.30  report  that  the  average 
composition  is  approximately  51  percent  quartz,  21  percent 
rock  fragments,  and  28  percent  feldspars,  highly  altered 
volcanic  components  with  diagenetic  cement  dominantly 
calcite.  In  uncemented  rocks  of  the  tidal  facies,  porosity 
averages  about  0.22  and  permeability  is  in  the  range  of  tens  to 
a  few  hundred  millidarcies,  with  a  mean  of  110-140  md 
(higher  for  flood-dominated  deposits). 

XMT  of  the  Sandstone 

A  sample  was  impregnated  with  an  epoxy  resin  under  vacuum 
and  then  cored  to  a  length  of  25  mm  and  diameter  of  5  mm.  A 
5  mm  long  section  of  the  core  was  imaged  at  33.07  keV 
energy  at  the  13-BMD  tomography  beamline,  operated  by 
GeoSoilEnviroCARS  at  the  Advanced  Photon  Source.  Image 
reconstruction  was  performed  using  algorithms  developed  by 
GSECARS  to  convert  CT  attenuation  data  to  3D  volumetric 
data  The  resultant  3D  gray-scale  image  has  a  voxel  resolution 
of  7.63  microns.  Segmentation,  the  process  of  converting  a 
gray-scale  image  to  a  1  -bit  or  binary  image  by  separating  the 
image  into  two  populations  based  on  gray-scale  values,  was 
performed  using  the  indicator  kriging  technique31  that  is  part 
of  the  3DMA  software  package.  In  this  case,  the  pore  voxels 
are  assigned  a  value  of  0  and  the  grain  voxels  are  assigned  a 
value  of  1.  Figure  1  is  an  image  of  the  binary  volume  file. 

While  the  porosity  of  the  resultant  binary  image  is 
relatively  close  to  that  of  the  bulk  sandstone  sample  (18.6% 
versus  22%),  the  combination  of  mineralology  (i.e,.  higher 
absorbing  elements)  and  X-ray  energy  created  a  relatively 
high  signal-to-noise  ratio.  This  less-than-ideal  image  quality 
made  it  difficult  to  completely  resolve  the  solid  and  void 
phases  and  to  remove  the  ring  artifacts.  Additional  imaging 
experiments  are  planned  that  should  produce  higher  quality 
datasets.  (See  concluding  comments.) 


A 


SPE  95887 


Figure  1.  Binary  volume  file  of  the  sandstone  obtained 
from  XMT 

Network  Construction 

Overview 

The  input  for  the  algorithm  is  a  binary  volume  file  describing 
the  porous  material  of  interest:  voxel  labels  indicate  whether 
the  voxel  is  contained  in  the  void  phase  or  solid  phase  of  the 
material.  The  algorilhim  operates  by  identifying  grains, 
searching  for  pores  based  on  the  grain  locations,  and  then 
creating  the  interconnected  network  using  a  novel  restricted- 
bum  algorithim.  The  significant  differences  between  this 
algorithm  and  other  network  generation  techniques  include  the 
following: 

1.  The  first  step  in  the  algorithm  is  the  identification  of  the 
granular  structure,  which  is  less  sensitive  to  image 
resolution  than  processes  such  as  skeletonization. 
Consequently,  networks  created  for  the  same  material  at 
two  different  image  resolutions  have  very  similar  structure 
and  properties. 

2.  The  algorithm  uses  the  granular  structure  as  a  template  to 
help  define  pores.  This  means  that  the  computationally 
expensive  search  for  pores  becomes  linked  to  the  number 
of  grains  in  the  image  rather  than  the  number  of  voxels, 
and  the  computational  penalty  for  increasing  image 
resolution  is  less  severe  than  with  other  algorithms. 

3.  Because  the  grain  locations  provide  the  framework  for 
determining  pore  structure,  the  algorithm  operates  without 
ever  computing  a  skeleton  of  the  pore  space.  By 
eliminating  this  step,  a  number  of  problematic  issues  such 
as  nonuniqueness,  dependence  on  image  resolution,  and 
the  formation  of  internal  loops  are  avoided. 

3.  A  key  intermediate  step  is  employed  between  the  original 
binary  map  and  the  final  network  structure,  which  is  the 
assignment  of  an  integer  label  to  every  voxel  (both  void 
and  solid  phase)  that  denotes  the  grain  or  the  pore  to  which 
that  voxel  belongs.  The  reason  for  emphasizing  this 
intermediate  step  is  that,  once  it  is  accomplished,  the  tasks 
of  obtaining  statistical  information  and  constructing  the 
network  are  straightforward  and  unambiguous. 

Grain  Reconstruction  Algorithm 

Step  1.  Dual-Phase  Burn .  A  simultaneous  grain-phase  and 
void-phase  bum  is  performed  (using  the  terminology  of 
Lindquist  et  al.,10).  Voxels  in  the  solid  phase  are  labeled  with 


positive  integers  denoting  the  bum  level  and  voxels  in  the  void 
phase  are  labeled  with  the  negative  of  the  burn  level.  This 
convention  is  not  necessary  if  the  burn  map  remains  coupled 
with  the  binary  material  array;  however,  the  use  of  opposing 
signs  allows  the  bum  numbers  to  be  written  over  the  initial 
material  array  without  losing  the  phase  information. 

Step  2.  Location  of  Extrema.  A  search  is  performed  to  find 
local  maxima  in  the  burn  assignments.  These  local  extrema  are 
islands  of  one  or  more  voxels  that  are  surrounded  by  voxels 
with  lower  burn  numbers.  In  the  simplest  form  of  the 
algorithm,  the  local  maxima  are  taken  to  be  the  grain  centers. 
However,  refining  the  grain  locations  using  optimization  leads 
to  better  results.  Simultaneously,  the  local  minima  can  be 
found  (which  according  to  the  sign  convention  are  in  the  void 
phase).  In  practice,  this  step  is  skipped  to  reduce  computation 
time  because  the  void-phase  extrema  are  not  used  to  define  the 
pore  structure.  However,  the  minima  are  found  and  reported 
for  completeness  in  this  paper. 

Step  3.  Tessellation  of  the  grain  centers .  A  periodic  Delaunay 
tessellation  is  performed  using  the  locations  of  the  grain 
centers.32  The  purpose  of  the  tessellation  is  to  identify  likely 
pore  locations  based  on  the  granular  structure.  The  role  of  the 
tessellation  is  somewhat  subtle,  and  the  following  points  are 
worth  clarifying: 

1 .  In  contrast  to  other  techniques  where  the  Delaunay 
tessellation  is  used  to  define  pore  structure, 9,1 8,1 9,25  it  does 
not  influence  the  structure  of  the  pores  in  the  current 
algorithm.  In  fact,  the  only  restriction  that  it  imposes  is  on 
the  total  number  of  pores:  the  algorithm  will  not  find  more 
pores  than  the  number  of  tetrahedrons  in  the  tessellation. 
This  limitation  should  not  be  of  practical  consequence 
since  the  tendency  of  the  Delaunay  tessellation  is  to 
identify  too  many  pores  (by  splitting  single  pores  into 
multiple  tetrahedrons;  see  reference  19).  Furthermore,  this 
restriction  can  be  relaxed  if  necessary  at  the  expense  of 
increased  computation  time  (see  item  3). 

2.  Use  of  the  Delaunay  tessellation  to  aid  in  the  identification 
of  pore  locations  provides  a  dramatic  advantage  as 
compared  to  using  more  traditional  erosion/dilation 
techniques  because  the  number  of  pores  that  are  located  for 
a  given  porous  medium  is  relatively  insensitive  to  voxel 
resolution.  This  issue  is  demonstrated  later  in  the 
validation  section. 

3.  Finally,  it  should  be  pointed  out  that  the  tessellation  is  a 
valuable  but  not  essential  part  of  the  algorithm.  An 
alternative  approach  would  be  to  perform  the  optimizations 
described  in  the  next  step  beginning  with  every  void-phase 
voxel.  However,  even  for  the  relatively  small  sandstone 
dataset  used  here,  this  approach  requires  solving 
-5, 000, 000  nonlinear  optimization  problems  (in  contrast  to 
-25,000  optimizations  when  the  tessellation  is  used  as  a 
template). 

Step  4 .  Locating  pores.  In  the  introduction  we  note  that  the 
division  of  void  space  into  pores  is  somewhat  arbitrary.  In  this 
work,  we  use  the  distance  function  d(x,y,z ),  whose  value  gives 
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the  minimum  distance  from  the  point  [x,y,z]  in  the  void  space 
to  a  grain  surface,33  and  then  define  a  pore  as  any  local 
maxima  in  d.  In  practice,  this  definition  corresponds  to  a 
common  definition  of  a  pore:  locations  of  maximum-diameter 
spheres  that  can  be  inscribed  into  the  void  space,  and  that  are 
constrained  from  movement  by  the  surrounding  solid  phase.34 
In  step  4,  these  maximum-diameter  inscribed  spheres  are 
located  by  performing  repeated  optimization  procedures  (to 
maximize  d)  using  the  tetrahedrons  as  seed  locations  to  start 
the  optimizations. 

The  optimizations  themselves  are  performed  using  a 
modified  Powell’s  method,35  which  is  a  direction  set  method 
effective  for  situations  where  gradients  in  the  objective 
function  cannot  be  calculated  directly.  In  essence,  the 
procedure  repeatedly  performs  one-dimensional  line 
minimizations.  Various  schemes  are  available  to  choose  the 
minimization  directions,  most  based  on  estimating  new 
conjugate  directions  from  the  history  of  the  optimization. 

As  each  extremum  is  found,  the  voxel  containing  its  xyj. 
location  is  marked  with  the  pore  number.  If  a  different  seed 
has  already  converged  to  this  same  voxel  (though  the 
coordinate  location  would  rarely  be  exactly  the  same),  a  new 
pore  is  not  added  to  the  list.  Hence,  Step  4  ends  having 
generated  a  list  of  N  pore  locations  and  inscribed  radii  (i.e., 
maximum  d  values),  along  with  N  void-phase  voxels  labeled 
with  the  corresponding  pore  number.  At  this  point,  the  pores 
are  no  longer  tied  to  their  seed  tetrahedrons,  and  there  is  no 
further  need  for  the  Delaunay  tessellation  in  the  algorithm. 

Step  5.  Pore  Merging.  A  viable  option  is  to  use  all  pore 
locations  identified  in  Step  4,  and  proceed  with  constructing 
the  network.  However,  in  real  cases,  many  pores  overlap  with 
one  another  by  a  significant  amount,  the  overlaps  being  caused 
by  one  of  two  reasons:  1 .  the  local  pore  geometry  causes  two 
independent  extrema  to  be  in  proximity;  2.  two  different  seeds 
have  led  to  the  same  extremum,  but  numerical  error  or 
optimization  tolerances  have  caused  the  computed  extrema 
locations  to  differ  by  at  least  one  voxel. 

In  either  of  the  above  cases,  there  is  good  reason  to  merge 
two  largely  overlapping  pores  into  one.  Various  merge  criteria 
can  be  devised.  In  the  current  algorithm,  pores  are  merged 
only  if  one  inscribed  sphere  encompasses  the  center  of  a 
neighboring  inscribed  sphere.  The  location  and  radius  of  the 
larger  inscribed  sphere  is  used  as  a  seed,  and  a  local 
optimization  is  performed  once  again  to  verify  the  location  of 
the  merged  pore.  With  the  adjusted  pore  locations  (due  to  the 
reoptimization  of  merged  pores),  we  have  found  it  advisable  to 
make  another  pass  to  check  whether  additional  pores  should 
be  merged,  and  indeed  to  continue  this  merge  — >  re¬ 
optimization  procedure  iteratively  until  no  more  pores  are 
merged.  Step  5  ends  with  the  same  information  as  Step  4  (pore 
locations,  radii,  and  corresponding  voxel  assignments),  but 
with  fewer  and  more  spatially  independent  pores. 

Step  6.  Grain  and  pore  assembly.  Step  6  is  the  key 
intermediate  step  mentioned  in  the  overview  above:  the 
assignment  of  all  voxels  in  the  image  to  one  of  the  grains  or 
pores  identified  in  Steps  2  or  5  respectively.  This  step  is 
performed  using  a  novel  restricted  burn  algorithm,  which  is 
described  in  more  detail  in  the  context  of  grain  reconstruction 


elsewhere.36  In  short,  is  assmbles  a  cluster  of  voxels  together 
that  are  tied  to  one  of  the  grain  or  pore  locations  respectively, 
in  accordance  with  the  local  geometry. 

A  summary  of  the  logic  for  this  restricted  bum  procedure 
(Step  6)  is  the  following: 

1 .  Set  the  minimum  bum  level  to  the  largest  (absolute)  value 
found  in  the  pore  space. 

2.  Loop  through  all  voxels  in  the  domain. 

3.  If  the  current  voxel  borders  a  voxel  already  assigned  to  a 
particular  pore  and  its  absolute  burn  level  is  greater  than  or 
equal  to  the  current  minimum  burn  level,  assign  it  to  the 
same  pore  as  the  assigned  neighbor. 

4.  If  any  new  voxels  were  assigned  during  the  last  pass,  go  to 

(2). 

5.  If  no  new  voxels  were  assigned  during  the  last  pass,  reduce 
the  minimum  bum  level  by  1;  go  to  (2). 

Step  7.  Pore  Morphology  and  Network  Construction.  Once 
Step  6  is  complete,  determining  morphologic  parameters  and 
constructing  the  network  is  a  straightforward  process.  The 
total  volume  of  each  pore  is  obtained  by  summing  the  volume 
of  all  voxels  assigned  to  that  particular  pore.  The  inscribed 
pore  radii  are  already  known  from  the  Step-4  and  Step-5 
computations. 

In  our  definition  of  network  structure,  pore  throats  have  no 
volume  but  rather  are  defined  by  the  faces  where  two  pore- 
elements  come  into  contact.  Hence,  the  pore-throats,  like  the 
pores,  have  rigorous  geometric  parameters  associated  with 
them.  The  unit  normal  for  the  pore-throat  interface  is  found  by 
averaging  the  orientation  of  all  voxel  interfaces  associated 
with  the  given  throat.  The  total  cross-sectional  area  of  the 
throat  is  then  found  by  summing  the  areas  of  voxel  faces  at 
this  interface  projected  onto  the  local  unit  normal.  (Use  of  the 
projected  area  prevents  overestimating  the  area  because  of  the 
staircase-like  convolutions  on  the  voxel  surface.)  The 
inscribed  area  of  the  throat  is  found  by  determining  the  largest 
inscribed  sphere  whose  center  is  located  on  the  throat 
interface. 

Grain  surface  area  is  also  assigned  as  a  throat  parameter 
since  it  affects  permeability  as  well  as  phenomena  such  as 
absorption  and  chemical  reactions  on  grain  surfaces.  The 
surface  area  is  assigned  by  estimating  the  surface  area  for  each 
surface  voxel,  and  then  assigning  that  element  of  surface  to  its 
closest  throat.  This  approach  ensures  that  the  surface  area  is 
conservative  (i.e.,  the  sum  of  all  pore  throat  surface  areas 
equals  the  total  surface  area  in  the  packing),  thus  providing  a 
good  theoretical  foundation  for  its  use  in  modeling. 

The  last  parameter  that  should  be  mentioned  is  the  pore- 
throat  hydraulic  conductance  of  each  throat,  which  is 
necessary  for  computing  permeability  and  performing 
dynamic  flow  simulations.  In  general,  the  conductance  is 
computed  using  some  combination  of  the  above-mentioned 
parameters,  via  generation  of  an  equivalent  capillary.  In  this 
work  we  use  the  same  approach  used  in  Thompson  and 
Fogler.18  Further  work  is  ongoing  to  evaluate  the  best  way  to 
compute  throat  conductivities  for  various  types  of  materials. 


The  network  itself  is  defined  by  mapping  out  the 
connectivity  of  each  pore,  which  is  defined  by  the  list  of 
neighboring  pores  that  share  voxel-voxel  contacts.  There  are 
no  limitations  on  the  structure  of  the  networks  thus  obtained 
(i.e.,  coordination  numbers,  etc.),  and  the  resulting  network 
files  have  exactly  the  same  format  as  network  files  from 
computer-generated  media.23  Once  the  network  is  constructed, 
plots  of  pore-size  distribution,  throat-size  distribution, 
coordination  number,  etc.,  can  be  made  directly  from  the  data 
in  the  network  file,  without  having  to  return  to  the  large 
voxelized  data  sets. 

Validation  (can  you  add  units  for  VR  etc.?) 

The  grain-reconstruction  algorithm  has  been  tested  extensively 
using  computer-generated  packings  of  spheres  and  cylinders 
as  well  as  XMT  data  from  unconsolidated  particles  and  sands. 
Separately,  the  network-generation  algorithm  has  been  tested 
using  a  series  of  computer-generated  packings  of  spheres.  This 
approach  is  valuable  for  three  reasons:  First,  the  pore  structure 
in  the  material  is  known  exactly;  hence,  the  validity  of  the 
network  generation  process  can  be  assessed  in  a  quantitative 
manner.  Second,  the  computer-generated  structures  can  be 
discretized  at  arbitrary  voxel  resolution,  which  provides  good 
benchmarks  for  the  accuracy  that  can  be  expected  with  data 
from  real  materials.  Third,  the  computer-generated  data  are 
free  from  noise  and  artifacts,  which  allows  validation  to  be 
focused  solely  on  the  network  construction  (rather  than 
imaging  and  segmentation  issues). 

Three  sphere  packings  were  used  for  validation:  a  cubic 
packing,  a  rhombohedral  packing,  and  a  random  packing.  The 
two  regular  packings  span  the  full  range  of  attainable 
porosities  for  monodisperse  spheres  and  the  pore  structure  is 
known  exactly.  The  random  packing  is  more  representative  of 
real  materials.  No  single  network  is  correct  for  the  random 
case.  However,  the  modified  Delaunay  Tessellation  (MDT) 


Table  II.  Notation  used  to  quantify  network  generation 
algorithms. _ 


Parameter 

Description 

Ars 

Cross  sectional  area  of  a  pore  throat  (ave) 

A, 

Surface  area  assigned  to  a  pore  throat  (ave) 

Dr. 

Grain  diameter  (ave) 

Dy 

Pore  diameter  (inscribed)  (ave) 

Dn 

Pore-throat  diameter  (inscribed)  (ave) 

H  EG 

Number  of  extrema  in  the  grain-phase  bum 

UEV 

Number  of  extrema  in  the  void-phase  bum 

Lht 

Length  of  pore  throat  (ave) 

MBG 

Maximum  bum  number  in  the  grain  phase _ 

MB  V 

Maximum  (absolute)  bum  number  in  the  void  phase 

Nr, 

Number  of  grains 

Nh 

Number  of  pores 

PLE 

Average  error  (%)  in  the  computed  pore  locations 

HP, 

Number  of  pores  initially  found  from  the  tetrahedron  seeds 

HPM 

Number  of  one-to-one  pore  matches  (reconstructed  packing  vs  original) 

HP,, 

Final  number  of  pores  after  the  iterative  merging/reoptimization  process 

n  Tct 

Number  of  tetrahedrons  used  as  seeds  for  pore  locations 

Vh,m„, 

Minimum  number  of  voxel  faces  required  for  a  pore  throat  to  be  formed 

H  Vox 

Number  of  voxels  in  the  data  set 

v„ 

Volume  of  a  pore  (ave) 

VPD 

Voxels  per  particle  diameter 

VR 

Voxel  resolution 

Z 

Pore  coordination  number 

£ 

Porosity 

structure  and  can  be  used  for  comparison  to  the  networks 
generated  from  the  voxelized  images. 

Table  2  explains  the  notation  used  for  reporting  the 
statistical  results.  Selected  quantitative  data  are  presented  in 
Table  3.  Note  that  for  the  three  sphere  packings,  the  top  row  of 
each  section  contains  parameters  obtained  from  the  MDT 
networks.  These  parameters  are  used  to  quantify  the  error  in 
the  voxel-based  parameters.  For  the  cubic  and  rhombohedral 
packings,  the  MDT  values  agree  with  theoretical  values  from  a 
unit-cell  analysis,  and  are  exact  within  the  numerical 
tolerances  set  in  the  MDT  algorithm.  As  mentioned  above, 
there  exist  no  “correct”  morphologic  values  for  the  random 
packing,  but  the  MDT  algorithm  is  viewed  as  an  excellent 
benchmark. 

Identification  of  Pores.  Column  #EV  is  the  number  of 
extrema  in  the  void-phase  burn.  This  parameter  is  the  most 
logical  choice  for  determining  pore  locations  if  the  bum 
information  were  to  be  used  directly.  However,  notice  the  very 


Table  III.  Validation  of  the  grain-based  network  generation  algorithm  using  computer-generated  sphere  packs. 


IMAGE  RESOLUTION 

BURN  PARAMETERS 

PORE  PARAMETERS 

PORE-THROAT  PARAMETERS 

VR 

VPD 

H  Vox 

MBG 

MBV 

HEG 

UEV 

H  Tel 

HP, 

HP,, 

L 

I7' 

Vr 

HPM 

PLE 

Z 

Drr 

A 

r.s 

Ax 

ave  | 

1  sd 

ave 

1  sd 

min 

|  max 

|  ave 

ave 

1  sd 

ave  | 

\sd 

ave  I 

sd 

Cubic  Packing 

MDT 

64 

0.7317 

0 

0.4764 

0 

6 

6 

6 

0.4142 

0 

0.5227 

0 

1.0472 

0 

0.100 

10.0 

40-40 >40 

3 

2 

64 

64 

458 

215 

64 

0.6996 

0.0003 

0.4480 

0 

64/64 

0.07% 

6 

6 

6 

0.3772 

0 

0.5046 

0 

0.9037 

0.0367 

0.096 

10.3 

41 -41 -41 

3 

3 

64 

64 

443 

201 

64 

0.7141 

0.0308 

0.4512 

0.0526 

64/64 

6.93% 

6 

13 

6.94 

0.3960 

0.0886 

0.4506 

0.1416 

0.7334 

0.2049 

0.090 

11.0 

44.44.44 

4 

2 

64 

64 

458 

218 

64 

0.7190 

0.0104 

0.4583 

0.0248 

64/M 

0.09% 

6 

6 

6 

0.4650 

0.0454 

0.5007 

0.0221 

0.8442 

0.0443 

0.080 

12.5 

50.50-50 

4 

3 

64 

64 

459 

139 

64 

0.7190 

0.0302 

0.4788 

0.0120 

64/64 

0.08% 

6 

6 

6 

0.4519 

0.0260 

0.5357 

0.0157 

0.9087 

0.0345 

0.050 

20.0 

80.80-80 

6 

4 

64 

64 

458 

270 

64 

0.7428 

0.0000 

0.4720 

0 

M/M 

0.06% 

6 

6 

6 

0.4078 

0 

0.5171 

0 

0.9247 

0.0404 

0.047 

21.3 

85.85*85 

7 

5 

64 

64 

418 

157 

64 

0.7189 

0.0119 

0.4742 

0.0082 

64/M 

2.37% 

6 

12 

6.66 

0.3856 

0.0978 

0.4766 

0.1399 

0.8594 

0.2069 

0.020 

50.0 

200-200.200 

14 

11 

64 

64 

458 

187 

64 

0.7366 

0.0001 

0.4740 

0 

M/M 

0.05% 

6 

6 

6 

0.4184 

0 

0.5166 

0 

0.9920 

0.0 1M 

0.010 

100.0 

400-400-400 

29 

21 

64 

64 

458 

194 

64 

0.7261 

0.0000 

0.4760 

0 

M/M 

0.05% 

6 

6 

6 

0.4198 

0 

0.5220 

0 

1 .0093 

0.0095 

Rhombohedral  Packing 

MDT 

192 

0.2878 

0.0892 

0.0612 

0.0498 

4 

8 

5J3 

0.1547 

0 

0.2266 

0 

0.3927 

0.0585 

0.100 

10.0 

40-34-32 

3 

2 

64 

4 

384 

277 

177 

0.3374 

0.0503 

0.0607 

0.0365 

167/192 

41.40% 

2 

20 

6.71 

0.2310 

0.0437 

0.1761 

0.0457 

0.2232 

0.0603 

0.060 

16.7 

66-57-54 

5 

2 

64 

160 

384 

261 

192 

0.2994 

0.0742 

0.0599 

0.0472 

192/192 

23.02% 

1 

20 

6.38 

0.1795 

0.0205 

0.1827 

0.0410 

0.2688 

0.0792 

0.040 

25.0 

100-86-81 

8 

4 

64 

88 

384 

271 

192 

0.2922 

0.0786 

0.0598 

0.0470 

192/192 

18.04% 

2 

19 

6.41 

0.1489 

0.0309 

0.1862 

0.0560 

0.2785 

0.0976 

0.020 

50.0 

200-173x163 

15 

6 

64 

192 

384 

320 

192 

0.2875 

0.0836 

0.0608 

0.0491 

192/192 

6.80% 

1 

26 

6.60 

0.1339 

0.0421 

0.1877 

0.0705 

0.2958 

0.1080 

0.015 

66.7 

226-230.217 

20 

9 

64 

251 

384 

295 

192 

0.2863 

0.0846 

0.0607 

0.0493 

192/192 

6.38% 

1 

30 

6.25 

0.1266 

0.0500 

0.1788 

0.0808 

0.2863 

0.1 161 

0.010 

100.0 

400-346-326 

29 

13 

64 

328 

384 

276 

192 

0.2901 

0.0999 

0.0605 

0.0493 

192/192 

4.14% 

1 

31 

6.10 

0.1301 

0.0505 

0.1885 

0.0785 

0.3057 

0.1094 

Random  Sphere  Packing 

MDT 

409 

.3591 

0.0921 

0.0784 

0.0782 

3 

12 

4.6 

0.2631 

0.0857 

0.3541 

0.0940 

0.3350 

0.1220 

0.100 

10.0 

43-43.43 

3 

3 

100 

63 

615 

489 

301 

.3975 

0.0946 

0.1063 

0.1069 

271/409 

40.47% 

2 

20 

6.71 

0.2900 

0.0767 

0.2620 

0.1214 

0.2471 

0.1296 

0.050 

20.0 

87-87.87 

6 

6 

100 

127 

619 

562 

346 

.3559 

0.0930 

0.0926 

0.1005 

311/409 

26.79% 

1 

20 

6.38 

0.2334 

0.0900 

0.2M9 

0.1351 

0.2563 

0.1386 

0.040 

25.0 

109-109-109 

8 

7 

100 

155 

619 

581 

369 

.3461 

0.0968 

0.0868 

0.1059 

310/409 

23.46% 

2 

19 

6.41 

0.2231 

0.0921 

0.2566 

0.1380 

0.2441 

0.1428 

0.025 

40.0 

175-175x175 

12 

11 

100 

278 

620 

599 

376 

.3402 

0.0942 

0.0853 

0.1088 

322/409 

23.40% 

1 

26 

6.60 

0.2083 

0.1015 

0.2484 

0.1479 

0.2392 

0.1448 

0.010 

100.0 

438-438-438 

29 

27 

100 

670 

619 

611 

401 

.3309 

0.0977 

0.0800 

0.1331 

334/409 

2 1 .44% 

1 

30 

6.25 

0.1937 

0.1100 

0.2240 

0.1583 

0.2424 

0.1518 

0.008 

125.0 

548-548-548 

37 

34 

100 

803 

619 

612 

407 

.3334 

0.0959 

0.0788 

0.1312 

350/409 

22.90% 

1 

31 

6.10 

0.1948 

0.1142 

0.2185 

0.1643 

0.2455 

0.1527 

algorithm19  is  a  farily  rigorous  method  for  extracting  the  pore 


strong  dependence  of  these  values  on  image  resolution  (at 
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least  for  the  non-cubic  structures).  For  the  rhombohedral 
packing,  the  number  of  void-phase  extrema  varies  from  4  to 
328  as  image  resolution  is  varied  over  an  order-of-magnitude 
(the  true  number  of  pores  is  192).  For  the  random  packing, 
#EV  varies  between  63  and  803  as  resolution  varies  from 
10  VPD  to  125  VPD  (the  best  estimates  for  number  of  pores  in 
this  packing  are  ~  400).  Clearly,  the  void-phase  bum  provides 
a  poor  indicator  of  pore  location,  which  is  related  to  the 
problems  with  skeletonization  that  were  mentioned  previously. 

In  the  current  algorithm,  the  grain  structure  is  used  as  a 
template  for  determining  pore  locations.  The  initial  seeds  for 
pore  locations  are  based  on  a  Delaunay  tessellation  of  the 
grain  locations.  Table  3  lists  the  number  of  Delaunay 
tetrahedrons  (see  column  UTet).  The  number  of  largest- 
inscribed-spheres  obtained  from  these  seeds  is  listed  in 
column  UP,.  Note  that  the  numbers  in  both  of  these  columns 
are  relatively  insensitive  to  image  resolution  since  they  are 
tied  to  the  grain  structure  rather  than  the  voxel  structure. 

Focusing  on  the  cubic  packing,  note  that  there  are 
significantly  fewer  inscribed  spheres  (i.e.,  pores)  found  than 
seeds.  The  reason  for  this  difference  is  that  to  assemble  a 
cubic  pore  requires  at  least  five  tetrahedrons.  If  the 
optimization  routine  were  exact,  all  five  seeds  would  converge 
on  the  same  central  pore  location.  In  reality,  more  than  one 
seed  will  usually  converge  on  the  same  central  voxel,  and  is 
listed  only  once.  The  fact  that  UP/  decreases  with  increasing 
resolution  is  somewhat  counterintuitive  (because  higher 
resolutions  provide  many  more  voxels  for  the  optimization  to 
land  in).  However,  the  reason  for  this  effect  is  the  increased 
accuracy  of  the  optimization  procedure  as  resolution  increases 
(which  is  a  consequence  of  more  precise  distance-to-surface 
calculations). 

The  more  important  value  is  the  final  number  of  pores 
after  merging  {UP^j).  For  the  cubic  packing,  this  value  is  exact 
for  all  resolutions.  For  the  rhombohedral  packing,  it  is  exact 
for  resolutions  of  16.7  VPD  and  above.  For  the  random 
packing,  an  exact  number  of  pores  can  never  be  defined 
unequivocally.  However,  the  values  obtained  from  the  current 
algorithm  and  the  MDT  algorithm  agree  reasonably  well. 
Results  from  the  random  packing  also  illustrate  the  efficiency 
of  the  grain-based  approach.  Examining  the  number  of 
tetrahedrons  (each  of  which  generates  a  seed  point  to  search 
for  a  pore)  shows  that  -  620  nonlinear  optimization 
procedures  were  performed  to  find  the  initial  pore  locations, 
independent  of  image  resolution.  In  contrast,  ~  32  million  of 
these  same  computations  would  have  been  required  for  the 
100  VPD  image  if  each  void-phase  voxel  were  tested  to  find 
its  maximum  inscribed  sphere. 

Accuracy  of  Pore  Locations.  The  pore  location  is  defined  as 
the  center  of  the  largest  sphere  that  can  be  inscribed  into  a 
given  void  space.  For  the  cubic  packing,  the  results  are 
essentially  exact  in  cases  where  an  integer  VPD  value  is  used. 
Even  in  the  non-integer- VPD  cases,  the  error  remains  low. 
The  rhombohedral  packing  is  a  much  more  rigorous  test  since 
the  pores  are  not  symmetric  with  respect  to  the  Cartesian 
voxel  grid  and  also  because  they  are  small.  The  error  in  pore 
locations  is  significant  at  low  resolutions  (where  pore 
dimensions  can  be  as  small  as  a  single  voxel),  but  it  appears  to 
decrease  monotonically  with  increasing  voxel  resolution. 


Table  IV.  Effect  of  constraining  pore-throat  connections 


PORK-THROAT  PARAMETERS 

1  Tmi„ 

Z 

Dpt 

Acs 

As 

min 

|  max 

|  ave 

ave  |  sd 

ave  | 

sd 

ave  |  sd 

MDT 

1 

12 

4.60 

0.2627  0.0846 

0.3548 

0.0963 

0.3369 

0.1233 

0 

1 

26 

6.60 

0.2083  0.1015 

0.2484 

0.1479  02392 

0.1448 

1 

1 

26 

6.60 

0.2083  0.1015 

0.2484 

0.1479  0.2392 

0.1448 

2 

1 

23 

6.08 

0.2215  0.0950 

0.2673 

0.1386  02607 

0.1475 

3 

1 

21 

5.87 

0.2269  0.0921 

0.2757 

0.1338 

0.2697 

0.1477 

4 

1 

20 

5.71 

0.2310  0.0898 

0.2820 

0.1299  0.2770 

0.1472 

5 

1 

20 

5.71 

0.2310  0.0898 

0.2820 

0.1299  02770 

0.1472 

10 

1 

19 

5.24 

0.2422  0.0842 

0.3019 

0.1164 

0.3018 

0.1414 

50 

1 

18 

4.67 

0.2523  0.0814 

0.3251 

0.1008 

0J386 

0.1255 

63 

0 

17 

4.57 

0.3113  0.0960 

0.2981 

0.1536  0.3108 

0.1792 

The  accuracy  for  the  random  packing  is  somewhat  harder 
to  assess,  mostly  because  the  number  of  pores  found  does  not 
agree  exactly  for  the  MDT  versus  voxel-based  networks. 
Column  UPM  in  Table  3  provides  statistics  for  the  fraction  of 
pores  with  a  one-to-one  match  with  the  MDT  network 
(defined  when  the  center  of  one  and  only  one  pore  from  the 
voxel  network  lies  inside  an  inscribed  pore  of  the  MDT 
network).  Also  shown  is  the  average  error  in  pore  location  for 
these  one-to-one  matches. 

Coordination  number.  Determining  the  pore  coordination 
number  (the  number  of  throats  emanating  from  a  pore  and 
connecting  to  other  pores)  remains  the  biggest  problem  for  the 
algorithm.  For  the  cubic  packing,  the  coordination  number  is 
exact  for  integer  VPD  values,  again  due  to  the  fact  that  pore 
boundaries  are  coincident  with  voxel  boundaries.  For 
resolutions  that  produce  non-integer  VPD  values,  the  average 
coordination  numbers  are  slightly  higher  and  the  maximum 
coordination  numbers  are  significantly  higher  than  they  should 
be.  The  problem  is  illustrated  in  Figure  2,  which  shows  two 
pores  assembled  from  a  cubic  packing  with  a  non-integer 
resolution  of  21.3  VPD.  Although  they  are  constructed 
correctly  (within  the  limitations  of  the  voxel  resolution),  the 
fact  that  the  pore  boundaries  are  not  coincident  with  voxel 
boundaries  creates  a  contact  between  two  voxel  faces  at  a 
position  that  should  be  a  comer-comer  contact  in  reality  (and 
therefore  should  not  register  as  a  throat  connection).  The 
numbers  from  the  non-integer  VPD  cases  indicate  that  this 
problem  occurs  less  often  than  once  per  pore  on  average 
(though  certain  pores  are  assigned  many  extra  neighbors,  as 
evidence  by  the  high  values  for  the  maximum  coordination 
numbers).  For  the  rhombohedral  packing,  the  problem  is 
similar,  though  no  resolution  generates  the  exact  coordination 
pattern  since  the  rhombohedral  pore  boundaries  never  lie 
along  voxel  boundaries. 


Figure  2.  Non-physical  pore  throat  connection  in  a  cubic 
packing  of  spheres. 
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For  random  pore  structures  the  problem  is  exacerbated. 
The  pore  geometries  are  more  irregular  (anomalies  are  not 
limited  to  comer-comer  contacts  as  with  the  cubic  packing) 
and,  like  the  rhombohedral  packing,  the  Cartesian  voxel 
system  does  not  allow  for  perfect  breaks  in  the  pore  geometry. 
Estimates  of  average  pore  coordination  number  fall  into  the 
mid-sixes  range.  These  values  are  probably  too  high,  again 
because  of  identification  of  voxel-voxel  connections  that 
should  not  register  as  pore-pore  connections  (a  fact  also 
evidenced  by  the  large  values  of  maximum  coordination 
number). 

The  obvious  solution  to  this  problem  is  to  set  a  minimum 
limit  for  the  number  of  voxel-voxel  connections  that 
constitutes  a  pore  throat.  For  the  cubic  packing,  requiring  pore 
throats  to  be  composed  of  three-or  more  voxel  faces  led  to 
perfect  coordination  numbers  at  all  non-integer  resolutions 
tested.  For  the  40  VPD  random  packing.  Figure  3  shows  a 
histogram  of  the  number  of  voxel  faces  found  in  a  throat 
connection.  The  histogram  has  a  broad  minimum  before  the 
population  begins  to  increase  at  around  70  voxel-faces/throat. 
Interestingly,  the  theoretical  minimum  for  throat  size  at  this 
resolution  is  64  voxel-faces/throat,  which  lends  credence  to 
the  histogram.  These  two  factors  suggest  that  connections 
made  at  less- than  ~64  voxel-faces/ throat  are  anomalies  and 
should  be  discarded.  Table  4  shows  the  results  for  various 
minimum- voxel  limits  applied  to  the  40  VPD  random  sphere 
packing.  For  a  conservative  minimum  of  50  voxel- 
faces/throat,  the  revised  throat  parameters  are  in  much  better 
agreement  with  the  MDT  values,  with  the  average  pore 
coordination  number  reduced  to  4.67.  A  more  aggressive  limit 
of  63  voxel-faces  per  throat  caused  the  values  to  deteriorate, 
indicating  that  real  throats  are  being  discarded  with  this  higher 
limit. 

Unfortunately,  it  is  difficult  to  generalize  an  effective  rule 
for  how  to  limit  interconnectivity.  Creating  a  Figure-3-type 
histogram  on  a  case-by-case  basis  is  a  sensible  and  fairly  easy 
approach  with  an  automated  algorithm.  However,  initial 
tomographic  data  from  real  materials  that  we  have  tested  do 
not  show  the  bimodal  distribution  found  in  Figure  3  and 
therefore  do  not  provide  a  strong  rationale  for  a  cutoff  value. 


Figure  3.  Histogram  of  the  number  of  pore-throat  voxels 
in  a  random  packing  of  spheres. 

We  suggest  the  following  approach.  For  the  purposes  of 
network  generation,  no  cutoff  should  be  used  because  extra 


throats  will  have  negligible  effect  on  most  transport  processes 
(due  to  their  very  small  size),  and  because  one  risks 
eliminating  what  the  tomography  has  identified  as  Teal’ 
connections.  For  the  purposes  of  statistical  analysis,  the  extra 
throats  do  cause  a  problem  because  they  affect  the  calculation 
of  coordination  number  and  average  throat  parameters.  In 
these  cases  it  is  worth  the  extra  effort  to  assess  the  distribution 
of  voxel  contacts,  and  provide  a  cutoff  if  a  case  can  be  made 
to  do  so.  A  fmal  point  is  that  the  use  of  geometric  averages 
should  improve  the  statistical  characterization  even  in  cases 
where  small  throats  are  mistakenly  identified. 

Results 

Two  different  network  reconstructions  have  been  performed 
on  the  volume  data  shown  in  Figure  1.  The  first  uses  the 
grain-reconstruction  algorithm.  For  comparison,  a  second  run 
was  performed  in  which  the  search  for  pores  was  performed 
directly  on  the  void-phase  voxels,  without  use  of  the  grain 
reconstruction  and  the  Delaunay  tessellation. 

All  runs  were  performed  on  IBM  Power5  8-way  p575 
machines  housed  at  LSU’s  high-performance  computing 
center. 

Morphologic  Parameters 

Figure  4  contains  an  image  of  the  grain  reconstruction,  which 
is  included  for  general  interest.  Colors  are  assigned  randomly 
to  grain  numbers  to  give  a  visual  indication  of  the 
discretization  into  individual  grains.  (The  missing  piece  at  the 
bottom  left  comer  is  an  artifact  caused  by  limitations  in  the 
graphics  software.) 

In  previous  tests  of  the  algorithm  on  tomography  images  of 
unconsolidated  sands,  two  problems  with  the  grain 
reconstruction  process  were  identified.  The  first  problem  is 
single  grains  being  “broken”  into  more  than  one  piece.  This 
problem  is  associated  with  noise  in  the  data  and/or  the 
misidentification  of  grain  centers  (which  in  turn  stems  from 
either  unusual  grain  shapes  or  limitations  in  the  voxel 
resolution  and  resulting  bum).  The  second  problem  is  the 
opposite  situation:  the  identification  of  a  single  large  grain  that 


Figure  4.  Grain  reconstruction  of  the  Figure-1  data  set. 

should  probably  be  broken  into  more  than  one  piece  (although 
this  decision  is  subjective  if  any  level  of  consolidation  is 
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Grain-Based 

Network 

Voxel-Based 

Network 

Na 

2,334 

A,  (/All) 

101 

nh 

10,768 

65.574 

DP(v m) 

34 

18 

DrAv m) 

31 

26 

LPr(pn\) 

77 

35 

Z 

3.21 

3.16 

e 

0.186 

0.187 

Kju  (mD) 

410 

430 

Kyy  (mD) 

530 

570 

K~  (mD) 

320 

360 

present).  Both  of  these  problems  appear  to  be  more 
pronounced  for  the  current  sandstone  data  as  compared  to  the 
cleaner  data  sets  available  for  the  unconsolidated  sands. 
Solutions  to  these  problems  have  been  implemented  in  the 
algorithm,  but  are  currently  undergoing  validation  tests  and  so 
were  not  used  in  the  current  work. 

Fortunately,  the  main  problem  with  the  current  image 
appears  to  be  the  division  of  single  grains  into  more  than  one 
piece,  which  will  not  affect  the  network  generation  process  in 
a  detrimental  way;  an  increased  number  of  grains  will  simply 
cause  the  algoriththm  to  search  for  more  pores.  However,  the 
final  pore  locations  are  determined  by  the  optimization 
procedure  and  reflect  local  geometry  regardless  of  whether  the 
local  grains  were  reconstructed  in  a  reasonable  way.  Hence, 
the  grain-based  algorithm  generates  a  physically 
representative  network  even  if  problems  occur  in  the  grain 
reconstruction  process. 

Table  5  contains  selected  parameters  generated  by  the 
algorithm  for  the  two  different  network  generation  techniques. 
Figure  5  shows  histograms  of  the  inscribed  pore-size, 
inscribed  pore-throat-size,  and  pore  coordination  number  for 
the  grain-based  network.  Figure  6  shows  network  structure 
from  the  same  central  region  of  the  sandstone,  but  for  the  two 
different  networks. 

The  most  striking  difference  between  the  two  networks  is 
the  number  of  pores,  with  the  grain-based  network  having 
many  fewer  pores  and  thus  longer  pore  throats  spanning  the 
gaps  between  the  pores.  This  is  a  consequence  of  the  logic  in 
the  grain-based  algorithm,  which  uses  the  Delaunay 
tessellation  of  the  grains  as  seed  locations  for  potential  pores. 
In  contrast,  the  voxel-based  algorithm  searches  for  pores  from 


inscribed  spheres).  Although  the  difference  in  pore  numbers  is 
extreme,  it  is  still  not  possible  to  declare  one  network  more 
correct  than  the  other  because  discretization  of  the  pore  space 
is  an  arbitrary  process  (within  reason  of  course).  The  key  is 
that  in  for  both  cases  the  network  is  a  one-to-one  mapping  of 
the  pore  structure  and  in  both  cases  the  pores  are  defined  using 
the  same  rigorous  definition:  a  local  maximum  in  the  distance- 
to-a-grain-surface  function.  The  most  effective  illustration  of 
these  points  is  to  examine  the  two  networks  side  by  side  as  in 
Figure  6  and  note  that,  despite  differences  in  the  details,  the 
trends  in  network  structure  are  the  same  (the  appearance  of 
large  pores  at  the  same  locations;  gaps  in  the  network  structure 
at  the  same  locations). 

The  other  notable  point  from  the  data  presented  here  is  the 
relatively  small  difference  between  the  pore-throat  diameters 
and  the  pore  diameters,  which  is  a  consequence  of  the  channel 
like  structure  in  consolidated  sands.  In  fact,  in  the  voxel-based 
network,  the  average  pore  size  is  smaller  than  the  average 
throat  size.  Although  this  is  counter  to  the  traditional  pore  and 
pore-throat  model,  it  is  simply  a  consequence  of  the  high 
density  of  pores  in  the  second  network.  This  high  pore  density 
causes  the  void  space  to  be  defined  as  strings  of  largely 
overlapping  pores  (rather  than  distint  pores  connected  by  long 
throats),  in  which  case  a  pore-throat  diameter  is  simply  the 
size  of  the  channel  at  the  point  where  two  inscribed  spheres 
overlap,  and  thus  is  not  necessarily  smaller  than  one  or  both  of 
the  adjacent  pores. 

Flow  Modeling 

The  two  networks  were  used  to  model  single-phase  low- 
Reynolds-number  flow  of  a  Newtonian  fluid.  The  flow 
modeling  was  performed  after  cropping  100  /jm  from  all  sides 
of  the  networks  so  as  to  avoid  edge  effects  that  are  caused  by 
the  boundary  of  the  volume  data. 

Pressure  gradients  were  imposed  in  one  of  the  coordinate 
directions  and  the  resulting  volumetric  flowrates  were 
computed.  Permeabilty  was  then  calculated  using  Darcy’s  law. 
The  ability  to  obtain  dimensional  permeability  values  is  a 
consequence  of  using  physically  representative  networks.  That 
is,  the  network  is  a  map  of  a  specific  region  of  the  sandstone, 
which  allows  for  the  computation  of  dimensional  volumetric 


Figure  5.  Morphology  of  the  pore  space  in  the  sandstone:  (a)  pore-size  distribution  (inscribed);  (b)  pore-coordination-number  distribution; 
(c)  pore-throat-size  distribution  (inscribed). 


each  voxel  location,  which  results  in  the  identification  of 
many  more  local  extrema  (i.e.,  locations  of  maximum 


flowrates,  cross-sectional  areas  for  flow,  etc.  Additionally, 
because  of  the  one-to-one  mapping,  the  model  has  no 


It 


Figure  6.  Comparison  of  a  section  of  the  grain-based  network  (left)  to  the  same  section  of  material  in  the  voxel-based  network  (right) 


adjustable  or  scaling  parameters.  (The  pore-throat  hydraulic 
conductivities  are  determined  using  the  formulas  developed 
for  sphere  packings  in  reference,18  and  thus  are  not  treated  as 
adjustable  parameters.) 

The  similarity  in  permeabilities  for  the  two  different 
networks  may  be  surprising  considering  the  different  network 
structures.  However,  this  fact  is  again  a  consequence  of  using 
physically  representative  networks:  in  the  network  with  fewer 
pores  and  pore  throats,  the  throats  are  shorter  and  therefore 
have  larger  hydraulic  conductivities.  Put  simply,  everything 
comes  out  in  the  wash. 

Figure  7  is  a  color  mapping  of  single  phase  flow  through 
the  grain-based  network  (cropped  as  described  above).  The 


Figure  7.  Color  representation  of  single-phase  flow  through  < 
physically  representative  network  model  of  the  sandstone. 

colors  are  mapped  according  to  volumetric  flowrate,  with 
brighter  colors  indicating  high  flowrates  and  the  blues 
indicating  low.  This  graphic  shows  clearly  the  heterogeneity  in 
the  flow  patterns  that  results  from  the  pore  morphology,  and 
demonstrates  the  rule  of  thumb  that  relatively  few  pores  carry 
a  large  fraction  of  the  fluid.  This  graphic  is  also  an  effective 
illustration  of  the  rationale  for  using  network  modeling  as 
opposed  to  simpler  effective  medium  or  bundle-of-tubes 


models.  Although  these  simpler  models  may  well  predict 
correct  permeabilities  given  the  proper  pore  statistics,  they  do 
not  capture  the  flow  heterogeneity  that  contributes  strongly  to 
processes  such  as  solute  dispersion. 

Concluding  Comments 

It  is  worth  reemphasizing  that  the  3D  image  data  used  in  this 
work  contains  more  noise  and  imaging  artifacts  than  we  are 
comfortable  with  in  the  context  of  performing  quantitative 
modeling.  Thus,  the  spirit  of  this  research  is  to  develop 
powerful  and  robust  methods  for  network  generation  from 
binary  volume  data,  and  to  understand  the  implications  of 
decisions  such  as  the  technique  used  for  pore  identification. 

Already,  x-ray  tomography  has  become  a  much  more 
accessible  technology  than  only  a  decade  ago,  and  there  is 
little  question  that  imaging  techniques  will  continue  to 
improve  rapidly.  This  underscores  the  importance  of 
developing  computational  tools  such  as  the  algorithm 
described  here. 
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